home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 8: LINUX Games / Linux Cubed Series 8 - LINUX Games.iso / games / video / fly8111-.000 / fly8111- / fly8 / ofplane.c < prev    next >
C/C++ Source or Header  |  1979-12-31  |  19KB  |  802 lines
  1. /* --------------------------------- ofplane.c ------------------------------ */
  2.  
  3. /* This is part of the flight simulator 'fly8'.
  4.  * Author: Eyal Lebedinsky (eyal@ise.canberra.edu.au).
  5. */
  6.  
  7. /* Dynamics of the floating point expermental model.
  8.  *
  9.  * This model will eventualy do a proper aerodynamic simulation. It will
  10.  * go as far as 6DOF (simplified) and will properly implement a linearized
  11.  * serodynamics model. This program is being modified regularly as I learn
  12.  * more about the subject and gather more data. If I get a stable version
  13.  * then it will be separated into it's own module so as not to constantly
  14.  * have a half working program. But not yet.
  15.  *
  16.  * Note that the external parameters (EP->*) and internal data (EX->*)
  17.  * are read as fixed point so some conversion factors are attached to
  18.  * each. Once converted these become SI units.
  19.  *
  20.  * ToDo:
  21.  * 1 use proper trans- and super-sonic calcs for Cd (wave drag).
  22.  * 2 add wing sweep to calcs?
  23.  *
  24.  *
  25.  * The dimentionless derivatives are
  26.  *
  27.  * CL        lift, calculated from wing data
  28.  * CD        drag, calculated from wing data
  29.  * Cy.dr    rudder sideforce
  30.  * Cy.beta    vx damping
  31.  *
  32.  * Cl.da    aileron effectiveness
  33.  * Cl.p        roll damping
  34.  * Cl.beta    dihedral effect
  35.  * Cl.dr    roll from rudder
  36.  * Cl.r        roll from yaw (always set to Cl/4.0)
  37.  *
  38.  * Cm0        total Cm at 0 alpha
  39.  * Cm.de    elevators effectiveness
  40.  * Cm.q        pitch damping
  41.  * Cm.alpha    alpha (stabilizer) induced pitch
  42.  * Cm.alphadot    alpha rate induced pitch
  43.  *
  44.  * Cn.dr    rudder effectiveness
  45.  * Cn.r        yaw damping
  46.  * Cn.beta    weathercock stability
  47.  * Cn.da    ailerons induced yaw
  48.  * Cn.p        roll induced yaw
  49.  *
  50.  *
  51.  * The aerodynamic coefficients are:
  52.  *
  53.  * t  = - Cd        *cos(beta)
  54.  * Cx =   Cl        *sin(alpha)    [lift componnent]
  55.  *      + t        *cos(alpha)    [drag componnent]
  56.  * Cz = - Cl        *cos(alpha)    [lift componnent]
  57.  *      + t        *sin(alpha)    [drag componnent]
  58.  * Cy = - Cd        *sin(beta)    [drag componnent]
  59.  *    + Cy.dr        *rudder
  60.  *      + Cy.beta    *beta
  61.  *
  62.  * Cl =   Cl.da        *ailerons
  63.  *      + Cl.beta    *beta
  64.  *    + Cl.dr        *rudder
  65.  *      + Cl.p        *p        *b/2V
  66.  *      + Cl.r        *r        *b/2V
  67.  *
  68.  * Cm =   Cm0
  69.  *    + Cm.de        *elevators
  70.  *      + Cm.q        *q        *c/V
  71.  *      + Cm.alpha    *alpha
  72.  *      + Cm.alphadot    *alphadot    *c/V
  73.  *      + Cn.beta    *beta
  74.  *    + Cn.da        *ailerons
  75.  *      + Cn.p        *p        *b/2V
  76.  *
  77.  * Cn =   Cn.dr        *rudder
  78.  *      + Cn.r        *r        *b/2V
  79.  *
  80.  * The aerodynamic forces are:
  81.  *
  82.  * X = Cx*q*S
  83.  * Y = Cy*q*S
  84.  * Z = Cz*q*S
  85.  *
  86.  * The aerodynamic moments are:
  87.  *
  88.  * L = Cl*q*S*b
  89.  * M = Cm*q*S*c
  90.  * N = Cn*q*S*b
  91.  
  92. */
  93.  
  94. #include "plane.h"
  95.  
  96. #include <math.h>
  97.  
  98.  
  99. #define FVONE    ((float)VONE)
  100. #define FFONE    ((float)FONE)
  101. #define F10FONE    (FFONE/10.0)
  102. #define FA2R    (C_PI/2.0/FFONE)    /* angles to radians */
  103. #define FR2D    (180.0/C_PI)        /* radians to degrees */
  104.  
  105. #define DOPITCH        (st.debug & DF_GPW)
  106. #define DOYPLANE    (st.debug & DF_GPY)
  107. #define FUNDAMENTAL    (st.debug & DF_GPZ)
  108.  
  109. static float fDD[] = {1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0, 1000000.0};
  110. #define CCfshow(n,t,v)    CCshow(p, (n), t, (long)(fDD[n]*(v)))
  111.  
  112. LOCAL_FUNC void NEAR
  113. fairdata (float height, float *rrho, float *rho, float *sos)
  114. {
  115.     int    irrho, irho, isos;
  116.  
  117.     airdata ((long)(height*VONE), 0, &irrho, &irho, &isos);
  118.     if (rrho)
  119.         *rrho = irrho/FFONE;
  120.     if (rho)
  121.         *rho  = irho/FFONE;
  122.     if (sos)
  123.         *sos  = isos/FVONE;
  124. }
  125.  
  126. /* Get the total weight, including fuel and load.
  127. */
  128. LOCAL_FUNC void NEAR
  129. calc_weight (OBJECT *p, float *weight)
  130. {
  131.     float    t, tt;
  132.  
  133.     tt = EP->weight + EX->fuel/100.0;
  134.     if ((t = EX->stores[WE_M61-1]) > 0)
  135.         tt += t * st.bodies[O_M61]->shape->weight/1000.0;
  136.     if ((t = EX->stores[WE_MK82-1]) > 0)
  137.         tt += t * st.bodies[O_MK82]->shape->weight/1000.0;
  138.     *weight = tt * 0.4536;                    /* lb->Kg */
  139. }
  140.  
  141. /* Get wing coefficients from foil parameters and flaps settings.
  142. */
  143. LOCAL_FUNC void NEAR
  144. fGetFoil (OBJECT *p, float Cl0, float max_cl, float min_cl, float flaps,
  145.     float leFlaps, float *a0, float *maxCl, float *minCl)
  146. {
  147.     float    t;
  148.  
  149.     *a0 = Cl0 + EP->FEff/FFONE * flaps;
  150.     t =  EP->FEffCl/FFONE/FA2R * flaps;
  151.     t += EP->lefEffCl/FFONE/FA2R * leFlaps;
  152.     *maxCl = max_cl + t;
  153.     *minCl = min_cl + t;
  154. }
  155.  
  156. /* Calculate Cl and Cdi.
  157. */    
  158. LOCAL_FUNC int NEAR
  159. fGetCl (OBJECT *p, float aoa, float a0, float pi_e_AR, float maxCl, float minCl,
  160.     float *Cl, float *Cdi, float *Fstall, float *Downwash)
  161. {
  162.     float    a, aoaPstall, aoaNstall, cl, fstall;
  163.     float    t;
  164.     int    stall;
  165.  
  166.     a = 1.0/(2.0*C_PI) + pi_e_AR;
  167.     aoa -= a0;                /* get absolute aoa */
  168.     aoaPstall = maxCl * a;
  169.     aoaNstall = minCl * a;
  170.     if (aoa > aoaPstall) {
  171.         cl = maxCl;
  172.         if (!(EX->flags & PF_NOSTALL)) {
  173.             t = (aoa - aoaPstall)*2;
  174.             if (t >= aoaPstall)
  175.                 fstall = 0.0;
  176.             else
  177.                 fstall = t / aoaPstall;
  178.         } else
  179.             fstall = 1.0;
  180.         cl *= fstall;
  181.         stall = 1;
  182.     } else if (aoa < aoaNstall) {
  183.         cl = minCl;
  184.         if (!(EX->flags & PF_NOSTALL)) {
  185.             t = (aoa - aoaNstall)*2;
  186.             if (t <= aoaNstall)
  187.                 fstall = 0.0;
  188.             else
  189.                 fstall = t / aoaNstall;
  190.         } else
  191.             fstall = 1.0;
  192.         cl *= fstall;
  193.         stall = 1;
  194.     } else {
  195.         fstall = 1.0;
  196.         cl = aoa / a;
  197.         stall = 0;
  198.     }
  199.     *Cl = cl;
  200.     *Cdi = pi_e_AR * cl * cl;
  201.     if (Fstall)
  202.         *Fstall = fstall;
  203.     if (Downwash)
  204.         *Downwash = cl * pi_e_AR;    /* or 2.0 * ? */
  205.     return (stall);
  206. }
  207.  
  208. /* We use SI units (newton, kilogram, meter, second, radian) but read in some
  209.  * data in other units (feet, pound, pound-force, degrees etc.).
  210. */
  211. extern void FAR
  212. dynamics_fplane (OBJECT *p, int action)
  213. {
  214.     POINTER    *ptr;
  215.     int    onground, stall, it;
  216.     float    t, weight, thrust;
  217.     float    rrho, rho, sos, sa, ca, sb, cb;
  218.     float    AR;        /* wing aspect ratio (1/AR) */
  219.     float    pi_e_AR;    /* 1/(e*pi*AR) */
  220.     float    aFlaps;        /* flaps */
  221.     float    leFlaps;    /* leading edge flaps */
  222.     float    alpha;        /* geometric angle of attack */
  223.     float    walpha;        /* wing angle of attack */
  224.     float    beta;        /* angle of sideslip */
  225.     float    alphadot;
  226.     float    downwash;    /* downwash angle */
  227.     float    Cl;
  228.     float    Cdi;
  229.     float    Cstall;        /* stall level, [1...0]=[none...full] */
  230.     float    Cd;
  231.     float    S;
  232.     float    b;
  233.     float    c;
  234.     float    V, h;
  235.     float    qS_V;
  236.     float    qS;
  237.     float    qS_V_b2;
  238.     float    qS_V_c;
  239.     float    lift, drag;
  240.     float    Crudder, Cailerons, Celevators;
  241.     float    uu, vv, ww;
  242.     float    pp, qq, rr;
  243.     float    XX, YY, ZZ;
  244.     float    LL, MM, NN;
  245.     float    Fx, Fy, Fz;
  246.     float    mx, my, mz;
  247.  
  248.     if (DOYPLANE) {
  249.         dynamics_yplane (p, action);
  250.         return;
  251.     }
  252.  
  253.     if (action)
  254.         return;
  255.  
  256.     if (dynamics_input (p))
  257.         return;
  258.  
  259.     ptr = p->pointer;
  260.  
  261.     onground = EX->flags & PF_ONGROUND;
  262.  
  263.     if (onground ? check_takeoff (p) : check_land (p)) {
  264.         p->flags |= F_HIT;
  265.         return;
  266.     }
  267.  
  268.     uu =  EX->v[Y] /FVONE;
  269.     vv =  EX->v[X] /FVONE;
  270.     ww = -EX->v[Z] /FVONE;
  271.  
  272.     pp =  p->da[Y] *FA2R*FVONE;
  273.     qq =  p->da[X] *FA2R*FVONE;
  274.     rr = -p->da[Z] *FA2R*FVONE;
  275.  
  276.     Celevators = -EP->MaxElevators*FA2R *
  277.                 (EX->elevators + EX->tElevators)/100.0;
  278.     Cailerons = -EP->MaxAilerons*FA2R * EX->ailerons/100.0;
  279.     Crudder = EP->MaxRudder*FA2R * (EX->rudder + EX->tRudder)/100.0;
  280.  
  281.     V = p->speed /FVONE;
  282.     h = p->R[Z] /FVONE;
  283.  
  284. CCfshow (1, "uu", uu);
  285. CCfshow (1, "vv", vv);
  286. CCfshow (1, "ww", ww);
  287. CCfshow (2, "pp", pp*FR2D);
  288. CCfshow (2, "qq", qq*FR2D);
  289. CCfshow (2, "rr", rr*FR2D);
  290. CCfshow (2, "elv", Celevators*FR2D);
  291. CCfshow (2, "ail", Cailerons*FR2D);
  292. CCfshow (2, "rud", Crudder*FR2D);
  293. CCfshow (1, "V", V);
  294. CCfshow (1, "h", h);
  295.  
  296. /* Automatically refuel when stationery on ground.
  297. */
  298.     if (onground && 0 == p->speed && 0 == EX->power)
  299.         supply (p, 0);
  300.  
  301.     fairdata (h, &rrho, &rho, &sos);
  302. CCfshow (4, "rrho", rrho);
  303. CCfshow (4, "rho", rho);
  304. CCfshow (1, "sos", sos);
  305.  
  306. /* These are common terms.
  307. */
  308.     S = EP->wing_area /FVONE;
  309.     b = EP->wing_span /FVONE;
  310.     c = EP->wing_cord /FVONE;
  311.  
  312.     qS_V    = rho * V / 2.0 * S;    /* q*S/V */
  313.     qS      = qS_V * V;        /* q*S */
  314.     qS_V_b2 = qS_V * b / 2.0;    /* q*S/V*b/2 */
  315.     qS_V_c  = qS_V * c;        /* q*S/V*c */
  316.  
  317. /* Engine thrust.
  318. */
  319.     f16engine (p, (short)(sos*FVONE));
  320.  
  321. /* Thrust is read as lbf/10.
  322. */
  323.     thrust = EX->thrust * 10.0 * 4.4482;            /* lbf->N */
  324. CCfshow (0, "thrust", thrust);
  325.  
  326. /* Current weight.
  327. */
  328.     calc_weight (p, &weight);
  329. CCfshow (0, "weight", weight);
  330.  
  331. /* Find our alpha, beta and friends.
  332. */
  333.     if (onground && fabs(uu) < 2.0)
  334.         alphadot = walpha = alpha = beta = 0.0;
  335.     else {
  336.         alpha = atan2 (ww, uu);
  337.         beta  = asin (vv/V);
  338.         walpha = alpha + EP->Aoffset*FA2R;    /* wing rigging */
  339.         alphadot = (-EX->a[Z]*uu - EX->a[Y]*ww)/FVONE
  340.                 / (uu*uu + ww*ww);
  341.     }
  342.     EX->aoa = (ANGLE)(walpha/FA2R);
  343.     sa = sin (alpha);
  344.     ca = cos (alpha);
  345.     sb = sin (beta);
  346.     cb = cos (beta);
  347. CCfshow (2, "alpha", walpha*FR2D);
  348. CCfshow (2, "beta", beta*FR2D);
  349. CCfshow (2, "alpha.", alphadot*FR2D);
  350.  
  351. /* Get flaps effect. We introduce automatic flaps in response to alpha.
  352. */
  353.     aFlaps = EX->flaps;
  354.     if (EX->flags & PF_AUTOFLAP) {
  355.         if (aFlaps) {
  356.             EX->leFlaps = EX->flaps;
  357.         } else {
  358.             if (walpha > EP->AFamin*FA2R) {
  359.                 aFlaps = (walpha-EP->AFamin*FA2R)
  360.                         * EP->AFrate/FFONE/FA2R;
  361.                 if (aFlaps > EP->AFmax*1.0)
  362.                     aFlaps = EP->AFmax*1.0;
  363.             }
  364.             if (walpha > EP->ALEFamin*FA2R) {
  365.                 EX->leFlaps = (int)((walpha-EP->ALEFamin*FA2R)
  366.                         * EP->ALEFrate/FFONE/FA2R);
  367.                 if (EX->leFlaps > 100)
  368.                     EX->leFlaps = 100;
  369.             } else
  370.                 EX->leFlaps = 0;
  371.         }
  372.     } else {
  373.         EX->leFlaps = 0;
  374.     }
  375.     aFlaps  = EP->MaxFlaps*FA2R * aFlaps/100.0;
  376.     leFlaps = EP->MaxLEFlaps*FA2R * EX->leFlaps/100.0;
  377.  
  378.     if (EX->flaps)
  379.         aFlaps = EP->MaxFlaps*FA2R * EX->flaps/100.0;
  380.     else if (walpha > EP->AFamin*FA2R && (EX->flags & PF_AUTOFLAP)) {
  381.         aFlaps = (walpha-EP->AFamin*FA2R) * EP->AFrate/FFONE;
  382.         aFlaps = EP->MaxFlaps*FA2R * aFlaps;
  383.         if (aFlaps > EP->AFmax*FA2R)
  384.             aFlaps = EP->AFmax*FA2R;
  385.     } else
  386.         aFlaps = 0.0;
  387.     leFlaps = 0.0;
  388.  
  389. /* Get some basic data.
  390. */
  391.     AR = EP->wing_area / (EP->wing_span * (float)EP->wing_span)*FVONE;
  392.                         /* 1/AR */
  393.     t = AR / C_PI;                /* 1/(pi*AR) */
  394.     pi_e_AR = t / (EP->efficiency_factor/100.0);
  395.                         /* 1/(e*pi*AR) */
  396. {
  397.     float    maxCl, minCl;
  398.     float    a0;
  399.  
  400.     fGetFoil (p, EP->Cl0*FA2R, EP->maxCl/F10FONE, EP->minCl/F10FONE, aFlaps,
  401.             leFlaps, &a0, &maxCl, &minCl);
  402.     stall = 0;
  403.     EX->flags &= ~PF_STALL;
  404.     if (fGetCl (p, walpha, a0, pi_e_AR, maxCl, minCl, &Cl, &Cdi, &Cstall,
  405.             &downwash)) {
  406.         EX->flags |= PF_STALL;
  407.         if (!(EX->flags & PF_NOSTALL))
  408.             stall = 1;
  409.     }
  410.     Cailerons *= Cstall;
  411. CCfshow (3, "Cl", Cl);
  412. CCfshow (3, "Cdi", Cdi);
  413. CCfshow (2, "e", downwash*FR2D);
  414. }
  415.  
  416. /* The ground effect formula is extracted from the graph in Smiths 'The
  417.  * Illustrated Guide To Aerodynamics' end of chapter 3.
  418. */
  419.     t = EP->wing_span /FVONE;
  420.     if (h < t) {                /* ground effect */
  421.         t = h / t;            /* h/b */
  422.         if (t < 0.1)
  423.             t *= 5.0;
  424.         else
  425.             t = 1.06 - 0.07 / t;
  426.         Cdi *= t;
  427. CCfshow (3, "Cdi", Cdi);
  428.     }
  429.  
  430.     Cd = EP->Cdp0/FFONE;
  431.     Cd += Cdi;
  432.     if (EX->airbrake)
  433.         Cd += EP->Cds/FFONE * EX->airbrake/100.0;
  434.     if (EX->equip & EQ_GEAR)
  435.         Cd += EP->Cdg/FFONE;
  436.     if (EX->stores[WE_MK82-1] > 0)
  437.         Cd += EX->stores[WE_MK82-1] * (EP->CdMK82/FFONE);
  438. CCfshow (3, "Cd", Cd);
  439.  
  440. /* We now calculate the new forces, based on the state of the plane.
  441.  * The given state is:
  442.  *  pp        Rotation around y (right wing down, roll)
  443.  *  qq        Rotation around x (nose up, pitch)
  444.  *  rr        Rotation around z (nose right, yaw)
  445.  *  Crudder    rudder left-turn angle
  446.  *  Cailerons    ailerons roll-right (cwise) angle
  447.  *  Celevators    elevators push-down angle
  448.  *  alpha    nose above wind angle
  449.  *  beta    nose left-of-wind angle
  450.  * The forces are:
  451.  *  XX        forward
  452.  *  YY        right
  453.  *  ZZ        down
  454.  * The moments are:
  455.  *  LL        around x (right wing down, roll)
  456.  *  MM        around y (nose up, pitch)
  457.  *  NN        around z (nose right, yaw)
  458.  *
  459.  * From each force we get acceleration [a = F/m] and then velocity [v += a*t]
  460.  * and position [x += (v + a*t/2)*t].
  461.  *
  462.  * From each moment we get angular acceleration [aa = M/I], then we get
  463.  * angular rate [da += aa*t] which we use to update the Euler angles.
  464. */
  465.  
  466. /* The auto-elevators reduces the elevators control sensitivity linearly
  467.  * with the velocity. It engages at the prescribed velocity. Note the fudge
  468.  * factor for the fundamental model.
  469. */
  470.     if (EX->flags & PF_AUTOELEVATOR) {
  471.         t = EP->AErate /FVONE;
  472.         if (FUNDAMENTAL) {
  473.             if (DOPITCH)
  474.                 t /= 1.8;    /* experimental */
  475.             else
  476.                 t /= 4.0;
  477.         } else
  478.             t /= 1.2;
  479.         if (V > t)
  480.             Celevators *= t/V;
  481.     }
  482.  
  483.     LL = MM = NN = 0.0;
  484.     XX = YY = ZZ = 0.0;
  485.  
  486.     
  487. /* Gravity.
  488. */
  489.     YY -= weight * GACC/FVONE * p->T[X][Z]/FFONE;
  490.     XX -= weight * GACC/FVONE * p->T[Y][Z]/FFONE;
  491.     ZZ += weight * GACC/FVONE * p->T[Z][Z]/FFONE;
  492.  
  493.  
  494. /* Engines thrust.
  495.  * If the engine is mounted at [x,z] (forward and above cg) at an angle of
  496.  * Ea above the x-y plane then we have:
  497.  * Eb = atan2 (z,-x)
  498.  * Er = sqrt (x*x+z*z)
  499. */
  500.     XX += thrust * cos (EP->Ea*FA2R);
  501.     ZZ -= thrust * sin (EP->Ea*FA2R);
  502.  
  503. /* Lift and drag.
  504. */
  505.     lift =  Cl * qS;
  506.     drag = -Cd * qS;
  507.  
  508.     Fy =  sb * drag;
  509.     t  =  cb * drag;
  510.     Fx =  sa * lift + ca * t;
  511.     Fz = -ca * lift + sa * t;
  512.  
  513.     YY += Fy;
  514.     XX += Fx;
  515.     ZZ += Fz;
  516. CCfshow (0, "lift", lift);
  517. CCfshow (0, "drag", drag);
  518. CCfshow (0, "X", Fx);
  519. CCfshow (0, "Y", Fy);
  520. CCfshow (0, "Z", Fz);
  521.  
  522. if (FUNDAMENTAL) {
  523. /*
  524.  **************************************************
  525.  *
  526.  * This section evaluates forces/moments from more
  527.  * fundamental aerodynamic principles.
  528.  *
  529.  **************************************************
  530. */
  531.  
  532. /* Lift induced pitching moment.
  533. */
  534.     my = -Fz*EP->ACy/FVONE - Fx*EP->ACz/FVONE;
  535.     MM += my;
  536. CCfshow (0, "Ml", my);
  537.  
  538. /* Wing pitching moment.
  539. */
  540.     my = qS * EP->Cm0w/FFONE;
  541.     MM += my * b;
  542. CCfshow (0, "Mw", my*b);
  543.  
  544. /* Tail (stabilizers and elevators).
  545.  * Here we actually simulate tailerons - the whole stabilizer is rotated
  546.  * as an elevator, not as tail flaps.
  547.  * Downash ignored.
  548. */
  549.     {
  550.     float    ClTail, CdTail, Tlift, Tdrag;
  551.     float    ta, sta, cta, dta;
  552.     float    TAR;
  553.     int    estall;
  554.  
  555.     ta = alpha - downwash;            /* actual alpha */
  556.     sta = sin (ta);
  557.     cta = cos (ta);
  558.     ta += EP->Toffset*FA2R;
  559.     ta += Celevators;            /* tailerons */
  560.  
  561.     if (uu)
  562.         ta += (dta = atan2 (-qq * EP->TACy/FVONE, uu));
  563.     else
  564.         dta = 0;
  565. CCfshow (2, "Tda", dta*FR2D);
  566. CCfshow (2, "Taoa", ta*FR2D);
  567.  
  568.     TAR = EP->tail_area / (EP->tail_span * (float)EP->tail_span)*FVONE;
  569.     t = pi_e_AR * TAR/AR;
  570.     estall = fGetCl (p, ta, 0.0, t, EP->TmaxCl/F10FONE, EP->TminCl/F10FONE,
  571.             &ClTail, &CdTail, 0, 0);
  572.     t = qS * EP->tail_area/EP->wing_area;
  573.     Tlift = ClTail * t;
  574.     Tdrag = -(CdTail + EP->Cdp0/FFONE) * t;
  575.     Fx =  sta * Tlift + cta * Tdrag;
  576.     Fz = -cta * Tlift + sta * Tdrag;
  577.     my = -Fz * EP->TACy/FVONE - Fx * EP->TACz/FVONE;
  578.  
  579.     MM += my;
  580.     XX += Fx;
  581.     ZZ += Fz;
  582. CCfshow (0, estall ? "Tstall" : "Tlift", Tlift);
  583. CCfshow (0, "Tdrag", Tdrag);
  584. CCfshow (0, "Xt", Fx);
  585. CCfshow (0, "Zt", Fz);
  586. CCfshow (0, "Mt", my);
  587.     }
  588.  
  589. /* Body drag. Quite a kludge.
  590. */
  591.     Fy  =    - qS * 10.0 * sb;
  592.     YY += Fy;
  593.     Fz  =    - qS * 0.5 * sa;
  594.     ZZ += Fz;
  595. CCfshow (0, "Dy", Fy);
  596. CCfshow (0, "Dz", Fz);
  597.  
  598. /* Engine thrust induced pitching moment.
  599. */
  600.     my = -thrust * sin ((EP->Ea + EP->Eb)*FA2R);
  601.     MM += my * EP->Er/FVONE;
  602. CCfshow (0, "Mp", my);
  603.  
  604. /* Pitching moment damping.
  605. */
  606. if (DOPITCH) {
  607.     my =      qS_V_c  * EP->Cmq/F10FONE   * qq;    /* damping */
  608.     MM +=  my * c;
  609. CCfshow (0, "Md", my*c);
  610. }
  611.     my =      qS_V_c  * EP->Cmalphadot/F10FONE * alphadot;
  612.     MM +=  my * c;
  613. CCfshow (0, "Ma.", my*c);
  614.  
  615. /* Rolling moment.
  616. */
  617. #define Clr    (Cl/4.0)
  618.     mx =    + qS      * EP->Clda/FFONE    * Cailerons
  619.         + qS_V_b2 * EP->Clp/FFONE     * pp    /* damping */
  620.         + qS_V_b2 * Clr               * rr    /* yaw induced */
  621.         + qS      * EP->Clbeta/FFONE  * beta;    /* hihedral */
  622.     LL += mx * b;
  623. #undef Clr
  624.  
  625.  
  626. /* Yawing moment.
  627. */
  628.  
  629. /* Rudder. We assume the foil is symmerical.
  630. */
  631.     {
  632.     float    ClRudd, CdRudd, Rlift, Rdrag;
  633.     float    ta, dta, maxCl, minCl, a0;
  634.     float    RAR;
  635.     int    rstall;
  636.  
  637.     ta = -beta;
  638.     if (uu)
  639.         ta += (dta = atan2 (-rr * EP->RACy/FVONE, uu));
  640.     else
  641.         dta = 0;
  642. CCfshow (2, "Rda", dta*FR2D);
  643. CCfshow (2, "Raoa", ta*FR2D);
  644.  
  645.     fGetFoil (p, 0.0, EP->RmaxCl/F10FONE, -EP->RmaxCl/F10FONE, Crudder,
  646.         0, &a0, &maxCl, &minCl);
  647.  
  648.     RAR = EP->rudd_area / (EP->rudd_span * (float)EP->rudd_span)*FVONE;
  649.     t = pi_e_AR * RAR/AR;
  650.     rstall = fGetCl (p, ta, a0, t, maxCl, minCl, &ClRudd, &CdRudd, 0, 0);
  651.     t = qS * EP->rudd_area/EP->wing_area;
  652.     Rlift = ClRudd * t;
  653.     Rdrag = -(CdRudd + EP->Cdp0/FFONE) * t;
  654.     Fx = sb * Rlift + cb * Rdrag;
  655.     Fy = cb * Rlift - sb * Rdrag;
  656.     mx =  Fy * EP->RACz/FVONE;
  657.     my = -Fx * EP->RACz/FVONE;
  658.     mz =  Fy * EP->RACy/FVONE;
  659.  
  660.     LL += mx;
  661.     MM += my;
  662.     NN += mz;
  663.  
  664.     XX += Fx;
  665.     YY += Fy;
  666. CCfshow (0, rstall ? "Rstall" : "Rlift", Rlift);
  667. CCfshow (0, "Rdrag", Rdrag);
  668. CCfshow (0, "Xr", Fx);
  669. CCfshow (0, "Yr", Fy);
  670. CCfshow (0, "Lr", mx);
  671. CCfshow (0, "Mr", my);
  672. CCfshow (0, "Nr", mz);
  673.     }
  674.  
  675. } else {
  676. /*
  677.  **************************************************
  678.  *
  679.  * This section evaluates forces/moments from the
  680.  * dimentionless derivatives.
  681.  *
  682.  **************************************************
  683. */
  684.     my =    + qS      * EP->Cmde/FFONE * Celevators;
  685.     MM +=  my * c;
  686. CCfshow (0, "Me", my*c);
  687.  
  688. /* Pitching moment damping.
  689. */
  690.     my =    + qS_V_c  * EP->Cmq/F10FONE   * qq;    /* damping */
  691.     MM +=  my * c;
  692. CCfshow (0, "Md", my*c);
  693.  
  694. /* Wing and tail stabilizer pitching moment.
  695. */
  696. if (DOPITCH) {
  697.     my =    + qS      * EP->Cm0/FFONE
  698.         + qS      * EP->Cmalpha/FFONE      * alpha
  699.         + qS_V_c  * EP->Cmalphadot/F10FONE * alphadot;
  700.     MM +=  my * c;
  701. CCfshow (0, "Mt", my*c);
  702. }
  703.  
  704. /* Rolling moment.
  705. */
  706. #define Clr    (Cl/4.0)
  707.     mx =    + qS      * EP->Clda/FFONE    * Cailerons
  708.         + qS      * EP->Cldr/FFONE    * Crudder
  709.         + qS_V_b2 * EP->Clp/FFONE     * pp    /* damping */
  710.         + qS_V_b2 * Clr               * rr    /* yaw induced */
  711.         + qS      * EP->Clbeta/FFONE  * beta;    /* hihedral */
  712.     LL += mx * b;
  713. #undef Clr
  714.  
  715.  
  716. /* Yawing moment.
  717. */
  718.     Fy  =      qS * (EP->Cydr/FFONE * Crudder + EP->Cybeta/F10FONE * beta);
  719.     YY += Fy;                    /* side force */
  720. CCfshow (0, "sf", Fy);
  721.  
  722.     mz =    + qS      * EP->Cndr/FFONE    * Crudder
  723.         + qS      * EP->Cnda/FFONE    * Cailerons
  724.         + qS_V_b2 * EP->Cnr/FFONE     * rr    /* damping */
  725.         + qS_V_b2 * EP->Cnp/FFONE     * pp    /* roll induced */
  726.         + qS      * EP->Cnbeta/FFONE  * beta;    /* sideslip induced */
  727.     NN += mz * b;
  728. CCfshow (0, "N", mz*b);
  729. }
  730.  
  731. EX->misc[5] = qS ? (int)(MM/(qS*c)*1000000.0) : 0;
  732. CCshow (p, 3, "Cm", EX->misc[5]);
  733.  
  734. {
  735.     MAT    II;
  736.     VECT    F, M;
  737.  
  738.     II[X][X] = (short)EP->Iyy;
  739.     II[Y][Y] = (short)EP->Ixx;
  740.     II[Z][Z] = (short)EP->Izz;
  741.     II[Y][Z] = (short)EP->Izx;
  742.  
  743.     F[X] =  (short)(YY/weight *FVONE);
  744.     F[Y] =  (short)(XX/weight *FVONE);
  745.     F[Z] = -(short)(ZZ/weight *FVONE);
  746.  
  747.     M[X] =  (short)(MM/weight /FVONE/FVONE/FA2R);
  748.     M[Y] =  (short)(LL/weight /FVONE/FVONE/FA2R);
  749.     M[Z] = -(short)(NN/weight /FVONE/FVONE/FA2R);
  750.  
  751.     LandGear (p, F, M);
  752.  
  753.     SixDOF (p, F, M, II);
  754. }
  755.  
  756.     Mobj (p);
  757.  
  758.     p->speed = ihypot3d (EX->v);
  759.     VMmul (p->V, EX->v, p->T);
  760.  
  761. #define MAX_SPEED    (1000*VONE)
  762.     if (p->speed > MAX_SPEED) {            /* temp */
  763.         it = fdiv (MAX_SPEED, p->speed);
  764.         EX->v[X] = fmul (it, EX->v[X]);
  765.         EX->v[Y] = fmul (it, EX->v[Y]);
  766.         EX->v[Z] = fmul (it, EX->v[Z]);
  767.         p->V[X]  = fmul (it, p->V[X]);
  768.         p->V[Y]  = fmul (it, p->V[Y]);
  769.         p->V[Z]  = fmul (it, p->V[Z]);
  770.         p->speed = ihypot3d (EX->v);
  771.     }
  772. #undef MAX_SPEED
  773.  
  774.     if (onground)
  775.         EX->flags &= ~PF_STALL;
  776.  
  777.     if (EX->Gforce > EP->max_lift || EX->Gforce < EP->min_lift)
  778.         EX->flags |= PF_GLIMIT;
  779.     else
  780.         EX->flags &= ~PF_GLIMIT;
  781.  
  782. /* Mach number.
  783. */
  784.     EX->mach = muldiv (p->speed, 1000, (int)(sos*FVONE));
  785.  
  786. /* pull up warning time.
  787. */
  788.     it = muldiv (10000, iabs(p->speed), 300*VONE);
  789.     it = muldiv (it, iabs(p->a[X]), D90);
  790.     if (it < 2000)
  791.         it = 2000;
  792.     EX->misc[8] = it;
  793. }
  794.  
  795. #undef FVONE
  796. #undef FFONE
  797. #undef FA2R
  798. #undef DOPITCH
  799. #undef DOYPLANE
  800. #undef FUNDAMENTAL
  801. #undef CCfshow
  802.